arXiv:1506.08765vl [cs.CV] 29Jun2015 


Spectral Motion Synchronization in SE(3) 


Federica Arrigoni, Andrea Fusiello Beatrice Rossi 

DIEGM - University of Udine AST Lab - STMicroelectronics 

Via delle Scienze, 208 - Udine (Italy) Via Olivetti, 2 - Agrate Brianza (Italy) 


Abstract 

This paper addresses the problem of motion synchro¬ 
nization (or averaging) and describes a simple, closed- 
form solution based on a spectral decomposition, which 
does not consider rotation and translation separately but 
works straight in SE(3), the manifold of rigid motions. Be¬ 
sides its theoretical interest, being the first closed form so¬ 
lution in SE(fi), experimental results show that it compares 
favourably with the state of the art both in terms of precision 
and speed. 

1. Introduction 

In this paper we address the motion synchronization 
(a.k.a. motion averaging or motion registration) problem 
in the Special Euclidean Group, SE{3), which consists 
in recovering n absolute motions, i.e. rigid 3D displace¬ 
ments expressed in an absolute (external) coordinate sys¬ 
tem, starting from a redundant set of relative (pairwise) 
motions. This problem appears in the context of structure- 
from-motion (SfM) - where the absolute motions represent 
orientations and positions of cameras capturing a 3D scene, 
and multiple point-set registration — that requires to find the 
rigid transformations that bring multiple 3D point-sets into 
alignment. 

In the literature on multiple point-set registration, the ori¬ 
gins of motion synchronization can be traced back to the 
frame space methods [23] that optimize the internal coher¬ 
ence of the network of rotations and translations applied to 
the local coordinate frames, as opposed to solutions that op¬ 
timize a cost function depending on the distance of corre¬ 
sponding points (e.g. [22, 6]). 

In the structure-from-motion literature, global methods, 
that first solve for the motion by optimizing the network of 
relative motions and leave the 3D structure recovery at the 
end, are fairly recent (e.g. [19]), although the origins of 
these approaches can be traced back to [10]. 

Almost all the techniques address the motion synchro¬ 
nization problem by breaking it up into rotation and trans¬ 
lation, and solving the two synchronization problems sepa¬ 


rately. 

For what regards rotation synchronization, a theoretical 
analysis of the problem is reported in [15]. The absolute 
rotations can be recovered by using the quaternion repre¬ 
sentation of SO{3), as done in [10], or by distributing the 
error over cycles in the graph of neighbouring views [23]. In 
[14] a cost based on the fi-norm is used to average relative 
rotations, where each absolute rotation is updated in turn 
using the Weiszfeld algorithm. Martinec in [19] casts the 
problem as the optimization of an objective function based 
on the ^ 2 -iionti of the compatibility error between relative 
estimates and unknown absolute orientations, and this ap¬ 
proach is extended in [2] where approximate solutions are 
computed either via spectral decomposition or semidefinite 
programming. The sum of unsquared deviations is proposed 
in [28] as a more robust self consistency error. Chatterjee in 
[8] exploits the Lie-group structure of rotations, and com¬ 
bines an ii averaging in the tangent space with an iteratively 
reweighted least squares (IRES) approach. In [4] the rota¬ 
tion synchronization problem is reformulated in terms of 
low-rank and sparse matrix decomposition. 

As for translation synchronization methods, a discrim¬ 
inating factor relevant to our analysis is whether they use 
only relative motion information, or, in addition, exploit 
image point correspondences (e.g. [19, 29]). Let us focus 
on the former, which are more similar to our SE{3) ap¬ 
proach. In [10] absolute positions are initialized as the least 
squares solution of a linear system of equations in the pair¬ 
wise directions and orientations, and they are then improved 
through IRLS. In ["] a fast spectral solution to translation 
synchronization is proposed by reformulating the problem 
in terms of graph embedding. The method presented in 
[2 1 ] first computes pairwise directions through a robust sub¬ 
space estimation and then derives absolute translations us¬ 
ing a semidefinite relaxation. In [ 18] a linear solution which 
minimizes a geometric error in camera triplets is presented, 
while in [20] relative translations are computed through an 
a-contrario trifocal tensor estimation, and then absolute po¬ 
sitions are recovered by using an £oo formulation. 

A different approach for motion synchronization is fol¬ 
lowed in [11] where rotations and translations are jointly 
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considered. This method exploits the Lie-group structure of 
SE{3) and uses an iterative scheme in which at each step 
the absolute motions are approximated by averaging rela¬ 
tive motions in tangent space. In [12] robustness is intro¬ 
duced through random sampling in the measurement graph. 
Originally proposed in the SfM framework, this technique 
was also applied to multiple point-set registration [13] and 
simultaneous localization and mapping (SLAM) [1], where 
additional information about data reliability is introduced in 
the form of covariance matrices. 

In this paper we propose a novel method for synchro¬ 
nizing relative motions in SE{3), based on computing the 
least four eigenvectors of a 4n x 4n matrix. Our approach 
is the hrst one that works in SE{3) and has a closed-form 
solution, being based on a spectral decomposition. It can be 
seen as the extension to SE{3) of the spectral synchroniza¬ 
tion proposed in [24] for SO{2) and generalized in [25, 2] 
to SO{3). 

The simple matrix formulation of our method leads im¬ 
mediately to a weighted formulation, in much the same way 
as [28] did for SO{N), that allows to embed it into an IRLS 
scheme in order to handle rogue measurements. 

Experimental results on synthetic and real data show that 
it compares favourably with the state of the art both in terms 
of accuracy and efficiency. 

2. Our Method 

The motion synchronization problem consists in recov¬ 
ering n absolute motions, i.e. rigid displacements in 
expressed in an absolute (external) coordinate system, start¬ 
ing from a redundant set of relative (pairwise) motions. 
Such relative information is usually cotTupted by a dif¬ 
fuse noise, in addition to sparse gross errors (outliers). Let 
S C {1, 2,.. ., n} X {1, 2,..., n} denote the set of available 
pairs, which can be viewed as the set of edges of an undi¬ 
rected hnite simple graph Q = (V, £"), where vertices in V 
correspond to absolute motions. In practical applications 
this graph is far from complete, due to the lack of overlap 
between some pairs of images/scans. However, there is a 
signihcant level of redundancy among relative motions in 
general datasets, which can be used to distribute the error 
over all the nodes, avoiding drift in the solution. 

Each motion can be viewed as an element of the Spe¬ 
cial Euclidean Group SE{3), which is the semi-direct prod¬ 
uct of the Special Orthogonal Group SO{3) with As 
a matrix group, SE{3) is a subgroup of the General Lin¬ 
ear Group GL(4), thus the inverse of a displacement and 
composition of displacements reduce to matrix operations. 
Accordingly, each absolute motion is described by a homo¬ 
geneous transformation 

M, = e SEi3) (1) 


where Ri G SO{3) and ti G represent the rotation 
and translation components of the i-th transformation. Sim¬ 
ilarly, each relative motion can be expressed as 

) e SEiS) (2) 

where RijGSO{3) and encode the transformation 

between frames i and j. The link between absolute and 
relative motions is encoded by the compatibility constraint 

My = M,M-i (3) 

which is equivalent to i?y = RiRj andty = — 
by considering separately the rotation and translation terms. 
Relative motions can be seen as measurements for the ratios 
of the unknown group elements. Einding group elements 
from noisy measurements of their ratios is also known as 
the synchronization problem [9, 24]. 

The remainder of this section is organized as follows. In 
Sec. 2.1 we describe properties that hold when all the rela¬ 
tive information is exact, necessary to dehne our technique. 
Then we derive our spectral solution to motion synchroniza¬ 
tion (Sec. 2.2). In Sec. 2.3 our method is embedded into an 
IRLS framework in order to handle outliers among relative 
motions. Einally, Sec. 2.4 briefly presents the extension of 
our method in SE{N). 

2.1. The Exact Case 

The absolute transformations can be recovered from (3) 
- up to a global motion - if we express it in a useful equiva¬ 
lent way that takes into account all the relative information 
at once. Eor simplicity of exposition, we first consider the 
case where all the pairwise motions are available. 

Let X G K 4 nx 4 n ygjjote the block-matrix containing the 
ideal (noise free) relative motions and let M G be 

the stack of the absolute motions, namely 



'Ml' 


( 

Mi2 

Mi„\ 

M = 

M2 

, = 

M21 

h . 
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M„_ 



Mn 2 

■ h ) 


where I 4 indicates the 4x4 identity matrix. If M~^ G 
]R4x4n jjjg concatenation of the inverse of absolute mo¬ 
tions, i.e. M~^ = [Mj“^ M 2 "^ ... M“^], then the 

compatibility constraint turns into X — MM~\ and hence 
rank(A) = 4. Note that here X is not symmetric posi¬ 
tive semideflnite, in contrast to the case of SO{3). Since 
M~^M = n/ 4 , we obtain 

XM = nM (5) 

which means that - in the absence of noise - the columns of 
M are 4 (independent) eigenvectors of X associated to the 
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eigenvalue n. Equation (5) is equivalent to 

(n/ 4 „ - X)M = 0. ( 6 ) 

Thus the columns of M are a basis for the 4-dimensional 
null-space of L = (n./ 4 „ — X). The matrix L resembles a 
block Laplacian, as it will we clarihed ahead. 

Conversely, any basis U for null(L) will not coincide 
with M in general, since it will not be composed of eu¬ 
clidean motions. Specihcally, it will not coincide with 
[0 0 0 1] in every fourth row. In order to recover M from 
U it is sufficient to choose a different basis for null(L) that 
satisfies such constraint, which can be found by taking a 
suitable linear combination of the columns of U. More 
precisely, let P G ]^nx 4 ra j^^trix such that 

PU G consists of the rows of U with indices multiple 
of four. The coefficient a,f3 G of the linear combina¬ 
tion are solution of 


in the the symmetric adjacency matrix A = [wij\. Ac¬ 
cordingly, the degree matrix D of the weighted graph is de¬ 
fined as Di^i = J2j St {ij)eS Equation ( 8 ) still holds 
with these dehnitions, thus our spectral method extends to 
weighted motion synchronization. 


2.2. Dealing with Noise 


We now consider the case where the pairwise motions 
are corrupted by noise, hence they do not satisfy equations 
(3) and ( 8 ) exactly. Thus the goal is to recover the absolute 
motions such that they are “maximally compatible” with the 
available relative information. In order to address this mo¬ 
tion synchronization problem, we consider an algebraic cost 
function that measures the residuals (in the Frobenius norm 
sense) of Equation ( 8 ), namely 


min 

MgSB(3)" 


LM 


2 


F 


(9) 


PUOL = 0, PUf3 = 1 (7) 

where the first equation has a three-dimensional solution 
space. Let dg be a basis for the null-space of PU. 

Thus the columns of M corresponding to rotations coincide 
(up to a permutation) with [Ua^,Ua^, U ag] and M is re¬ 
covered as M = [/[«!, OL^, CKg, (3]. 

Note that this post-processing on the eigenvectors is not 
required for the spectral method in SOii), since any or¬ 
thogonal basis for the null-space of L coincides (up to a 
permutation) with the stack of the absolute rotations. 

We now consider the case of missing data, in which the 
graph Q is not complete. In this situation missing pairwise 
motions correspond to zero blocks in X. Let A G 
be the adjacency matrix of Q and let D G be the 

degree matrix of 0 , i.e. the diagonal matrix that contains 
the degree of node i in its entry 17^ ^. It can be seen that 
Eq. ( 6 ) generalizes to 

((L> - A) 14x4) o X)M = 0 ( 8 ) 

where 0 denotes the Kronecker product and o denotes the 
Hadamard product. The matrix {D — A) is the Laplacian 
matrix of the graph Q, which gets “inflated” to a 4 x 4- 
block structure by the Kronecker product with l 4 x 4 (a ma¬ 
trix hlled by ones), and then is multiplied entry-wise with 
X. Thus the columns of M are a basis for the 4-dimensional 
null-space of L = {{D — A) 0 I 4 X 4 ) o X. 

If Q is complete then D = (n—1) J„ and A = l„xn~^n. 
hence the matrix L reduces to the previous one. 

Weighted graph. 

In some applications we are given non-negative weights Wij 
that reflect the reliability of the pairwise measurements. In 
other words, ^ is a weighted graph with real weights, stored 


with the additional constraint |jm 4 |l^ = c in order to fix 
the global scale. Here 1114 denotes the fourth column of 
M, X denotes a noisy version of the ideal matrix X, which 
contains the measured relative motions G SE(3), and 
L = ((D — A) 0 14 x 4 ) o X. Hereafter we will consistently 
use the hat accent to denote noisy measurements. Such a 
problem is difficult to solve since the feasible set SE{3)'^ = 
SE{3) X • • • X SE{3) is non-convex. 

In order to make the computation tractable, we do not 
solve Problem (9) directly, but we proceed as follows. First, 
we look for an orthogonal basis for the (approximated) 4- 
dimensional null-space of L, by solving the following opti¬ 
mization problem 


min 

C/T U=nli 



2 

F ' 


( 10 ) 


In other words, we solve the homogeneous system of equa¬ 
tions LU = 0 in the least-squares sense, where the solution 
space is known to have approximately dimension 4. 

Then, we hnd an estimate for M within this space by 
forcing the solution to coincide with [0 0 0 1 ] in every 
fourth row. Finally, we project in iS'0(3) all the 3x3 blocks 
corresponding to rotations by using Singular Value Decom¬ 
position (SVD). 

Proposition 1. Problem (10) admits a closed-form solution, 
which is given by the 4 eigenvectors of LAL associated to 
the 4 smallest eigenvalues. 

Proof. We hrst observe that Problem (10) coincides with 
min ti {U~^ (L^ L)U). (11) 

C/T U=nl4. 

Let P be the unconstrained cost function corresponding to 
this problem, namely 

P{U) = {L'^L)U) + tr{A{U'^U - n/ 4 )) (12) 
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where A € is a symmetric matrix of unknown La¬ 

grange multipliers. Setting to zero the partial derivatives of 

with respect to U we obtain 

BT ^ ^ 

— = 2(L^L)U -f 217yl = 0 ^ {L^L)U = -UA. (13) 
aU 

Let Ui be any four eigenvectors of LAL (normalized so that 
||ui|| = ^Jn) and let be the corresponding eigenvalues. 
Then U = [ui|u 2 |u 3 |u 4 ] satisfies both (13) and the con¬ 
straint U'^U = nii, with A = — diag(Ai, A 2 , A 3 , A 4 ) (in¬ 
deed LAL admits an orthonormal basis of real eigenvectors 
since it is symmetric). In other words, any quadruple of 
eigenvectors is a stationary point for the objective function 
F. The minimum is attained in (11) if u, are the 4 least 
eigenvectors of LAL. □ 

Proposition 1 guarantees that the solution to problem 
(10) is given by the 4 least eigenvectors of LAL, which co¬ 
incide with the 4 least right singular vectors in the Singular 
Value Decomposition (SVD) of L. Such a solution repre¬ 
sents the best (in the Frobenius norm sense) 4-dimensional 
approximation for null(L). Within such a space, we find 
the solution that is closest to have every fourth row equal to 
[0 0 0 1] by solving system (7) in the least-squares sense. 
Then, such a solution is projected onto SE{3A - as in [5] - 
by forcing every fourth row to [0 0 0 1 ] and projecting 3x3 
rotation blocks onto S'0(3) through SVD. 

This technique has the advantage of being extremely fast, 
as motion synchronization is cast to eigenvalue decomposi¬ 
tion of a 4n X 4n matrix. Moreover, in practical application 
the measurement graph Q is sparse, thus employing sparse 
eigen-solvers (such as Matlab eigs) increases its effi¬ 
ciency. From the computational complexity point of view, 
the Lanczos method (implemented by eigs) is “nearly lin¬ 
ear” since every iteration is linear in n, if the matrix is 
sparse, but the number of iterations is not constant. 

2.3. Dealing with Outliers 

The fact that our spectral method copes easily with 
weights on individual relative motions allows a straightfor¬ 
ward extension to gain resilience to rogue input measures 
via Iteratively Reweighted Least Squares (IRLS). 

First, we solve (10) to obtain an estimate for M with 
given weights' as explained in the previous section, then 
we update the weights using the current estimate of absolute 
motions, and these steps are iterated until convergence. In 
our experiments we used the Cauchy weight function [16] 


* The initial weights are all 1 by default, but they can be initialized from 
any reliability information coming from the relative motion estimation pro¬ 
cedure. 


where rij = \\Mij — MiM~^\\F- The tuning constant c 
have been chosen, as customary, based on the median ab¬ 
solute deviation (MAD): c = 1.482 0med(|r — med(r)|), 
where med() is the median operator, r is the vectorization 
of the residuals , and 9 = 2. 

2.4. Generalization to SE(N) 

In this paper we focus on SE{3) because this group 
arises in several applications. However, it is straightforward 
to see that our analysis and the derived spectral method ap¬ 
ply equally well to any dimension. 

Suppose that we are given a redundant number of pair¬ 
wise ratios Mij G SE{N), and we want to estimate the 
associated group elements Mi G SE{N), which repre¬ 
sent rigid displacements in If the graph is com¬ 

plete then - in the absence of noise - the block-matrix 
X G K(''v+i)"x(Ar+i)n N + 1, and the columns 

of M are V -b 1 eigenvectors of X with eigenvalue n. If 
the graph is not complete then Equation ( 8 ) still hold, and 
hence the columns of M form a basis for the {N + 1)- 
dimensional null-space of L. Thus we can generalize our 
spectral method to synchronize elements of SE{N), by 
computing the -b 1 least eigenvectors of L. 

3. Experiments 

In this section we evaluate our spectral method - hence¬ 
forth called EIG-SE(3) - on both simulated and real data 
in terms of accuracy, execution cost and robustness to out¬ 
liers. We compare EIG-SE(3) to several techniques from the 
state-of-the-art. All the experiments are performed in Mat- 
lab on a MacBook Air with i5 dual-core @1.3 GHz. In 
order to compare estimated and ground-truth absolute mo¬ 
tions, we find the optimal transformation that aligns them 
by applying single averaging [14] for the rotation term and 
least-squares for the scale and translation. We use the angu¬ 
lar distance and Euclidean norm to measure the accuracy of 
absolute rotations and translations respectively. 

3.1. Simulated Data 

In these experiments we consider n absolute motions in 
which rotations are sampled from random Euler angles and 
translation components follow a standard Gaussian distri¬ 
bution. The level of sparsity of the measurement graph 
is defined through the average degree d of nodes. The 
available pairwise motions are corrupted by a multiplica¬ 
tive noise, where the rotation component has axis uniformly 
distributed over the unit sphere and angle following a Gaus¬ 
sian distribution with zero mean and standard deviation 
(Jr G [1°, 10°], and the translation components are sampled 
from a Gaussian distribution with zero mean and standard 
deviation jr G [0.01,0.1]. In this way we perturb both 
direction and magnitude of pairwise translations. All the 
results are averaged over 50 trials. 
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Figure 1. Mean angular errors (degrees) on the absolute rotations with d = 5 (left) and d = 30 (right). The value of ctt is meaningless. 



Figure 2. Mean errors on the absolute translations with d = 5 (left) and d = 30 (right). 


We evaluate the effect of noise on rotations and transla¬ 
tions both separately and together, by considering n = 100 
absolute motions, in the cases d = 5 and d = 30, which 
correspond to about 95% and 70% of missing pairs, respec¬ 
tively. Higher values of d correspond to better conditioned 
problems, with the same qualitative behaviour as d = 30. 
Please note that in the real cases reported in Tab. 3, the per¬ 
centage of missing pairs ranges from 30% to 90%. 

Rotation synchronization. As for rotations, besides 
Govindu-SE(3) [11], we consider general synchronization 
techniques such as the Weiszfeld algorithm [14], spectral 
relaxation [2] (EIG), semidehnite programming [2] (SDP), 
the Ll-IRLS algorithm [8], and the R-GoDec algorithm 
[4]. Methods based on quaternions (such as [10]) have been 
already proved inferior to the other methods in [19]. The 
code of Ll-IRLS is available on-line, while in the other 
cases we used our implementation. 

Eigure 1 reports the mean angular errors on the absolute 
rotations as a function of an, obtained by running the rota¬ 
tion synchronization techniques mentioned above. The best 
accuracy is obtained by EIG-SE(3) together with EIG, SDP 
and Govindu-SE(3). On the contrary, the robust approaches 
R-GoDec, Ll-IRLS and Weiszfeld yield worse results, to 
different extents, because they inherently trade robustness 
for statistical efficiency. 

The noise on relative translations does not have any influ¬ 


ence on absolute rotations, hence the value of ar is mean¬ 
ingless in this experiment. 

Translation synchronization. As for translations, we 
consider only methods working in frame space, i.e. not 
requiring point correspondences, such as SDR [21], the 
graph-embedding approach by Brand et al. [i'] and the 
works of Govindu [10, 11]. Among these methods, only 
EIG-SE(3) and Govindu-SE(3) [11] are influenced by the 
noise on the translation norms, for they work in SE{3), 
while this does not influence the remaining algorithms, 
which take as input relative translation directions. The code 
of SDR is available on-line, while in the other cases we 
use our implementation. In this simulation we do not per¬ 
turb the relative rotations {an = 0), thus all the methods 
are given ground-truth relative/absolute rotations. Noise on 
rotational component influences also the translation errors 
with results qualitatively similar to those reported here. 

Eigure 2 shows the mean errors on the absolute transla¬ 
tions as a function of ar (units are commensurate with the 
simulated data), obtained by running the techniques men¬ 
tioned above. Both EIG-SE(3) and Govindu-SE(3) outper¬ 
form all the analysed methods in terms of accuracy. 

When the measurement graph is extremely sparse (d = 
5) the methods by Govindu [10] and Brand et al. [7] yield 
larger errors than usual; by inspecting the solution it is 
found that this corresponds to wrong solutions concentrated 
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Figure 3. Mean errors on the absolute translations with noise on both rotations and translations, with d = 5 (left) and d = 30 (right). 


around a few locations. This can be visualized in Fig. 4, 
which shows ground-truth and estimated positions (after 
alignment) for a single trial when ctt = 0.08. Such a 
behaviour - which is called “clustering phenomenon” - is 
analysed in [21] where the cause has been traced back to a 
lack of constraints on the location distances. For this reason 
the authors of [2 1 ] introduce ad-hoc constraints in the mini¬ 
mization problem, forcing the differences between locations 
to be “sufficiently” large. On the contrary, EIG-SE(3) and 
Govindu-SE(3), by working in SE{3), implicitly enforce 
such constraints as they take in input the relative transla¬ 
tions with their norm. 


nization, by varying the number of absolute motions from 
n = 100 to n = 1000, all the others parameters being fixed. 
More precisely, we choose the values d = 10, cry = 0.05 
and an = 5° to define sparsity and noise. Eigure 5 reports 
the running times of the analysed algorithms as a function 
of the number of nodes in the measurements graph, showing 
that EIG-SE(3) is remarkably faster than Govindu-SE(3) 
and SDR. Indeed SDR solves a semidefinite programming 
problem and Govindu-SE(3) uses an iterative approach in 
which absolute motions are updated by performing multi¬ 
ple averaging in the tangent space; both these operations are 
more expensive than computing the least four eigenvectors 
of a dn X 4n matrix. 


X Ground-truth 
O EIG-SE{3) 

+ Govindu-SE(3) 
* SOR 
V Brand at a 
Govindu 


— * 


Figure 4. Clustering phenomenon. By enlarging the figure, the 
reader will distinguish one cluster of black triangles (v) and one of 
cyan squares (□) near the origin, which correspond to the locations 
obtained by Brand et al. and Govindu respectively. 



n n 


Figure 5. Execution times (seconds) of motion synchronization. 
A magnification is shown on the right to appreciate the timing of 
EIG-SE(3). 


Motion synchronization. In this experiment we consider 
all the pipelines that cope with both rotation and transla¬ 
tion synchronization, and work in frame space, namely SDR 
[21] and Govindu-SE(3) [11]. SDR has an original transla¬ 
tion stage while rotations are computed by iterating the EIG 
method. Our approach and Govindu-SE(3) recover both ro¬ 
tations and translations at the same time. 

Figure 3 reports the mean errors on the absolute trans¬ 
lations obtained after perturbing both relative rotations and 
translations. All the methods return good estimates, which 
are further improved by increasing edge connectivity, and 
EIG-SE(3) together with Govindu-SE(3) achieves the low¬ 
est errors. 

We also analysed the execution time of motion synchro- 


The rundown of these experiments is that, EIG-SE(3) 
achieves the same optimal accuracy of its closest com¬ 
petitor [11] in considerably less time. 

OutUers influence. In this experiment we study the re¬ 
silience to outliers of EIG-SE(3) with IRES. We consider 
n = 100 absolute motions sampled as before and we fix 
d = 30 to define sparsity. Since we are interested in 
analysing exact recovery in the presence of outliers, noise 
is not introduced in this simulation. The fraction of wrong 
relative motions - randomly generated - varies from 10% 
to 50%. Figure 6 reports the mean errors obtained by EIG- 
SE(3) and its IRES modification: the empirical breakdown 
point of EIG-SE(3) H- IRES is about 45%. 
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tamination. 

3.2. Real Data 

We apply EIG-SE(3) with IRES to the structure from 
motion problem, considering both the EPEE benchmark 
[27] and unstructured, large-scale image sequences from 
[29]. The latter are available on-line together with the rela¬ 
tive motions, while for the EPEE benchmark we computed 
them following a standard approach based on the essential 
matrix with a final bundle adjustment (BA) refinement of 
camera pairs. 

Owing to the depth-speed ambiguity, the magnitude of 
relative translations (also referred to as epipolar scales) are 
undefined. Therefore, the input relative motions do not fully 
specify elements of SE{3), and the unknown scales have to 
be computed. 

A straightforward approach (suggested in [11]) consists 
in iteratively updating these epipolar scales, i.e. during each 
iteration the scale of the translation of is set equal to 
that of where and Mj are the current estimates 

of camera motions. The starting scales are all equal to 1 and 
the procedure is iterated until convergence. In our imple¬ 
mentation this is combined with IRES in the same loop: in 
one step we update the IRES weights and in the next step 
we update the epipolar scales. 

A different approach is proposed in [3], where a two- 
stage method is developed for computing the epipolar scales 
based on the knowledge of two-view geometries only. Eirst, 
a Minimum Cycle Basis (MCB) for the measurement graph 
is extracted by using Horton’s Algorithm [17], then all the 
scales are recovered simultaneously by solving a homoge¬ 
neous linear system. This approach is based on the obser¬ 
vation that the compatibility constraints associated to these 
cycles can be seen as equations in the unknown scales. In 
this way all the epipolar scales are computed before per¬ 
forming motion synchronization. 

However, computing the epipolar scales is not part of the 
synchronization task, strictly speaking. As a matter of fact, 
this indeterminacy is an idiosyncrasy of the structure from 
motion problem, which is not shared, e.g., by the multiple 
point-set registration problem, where the relative motions 
are fully specified. Eor this reason we are agnostic about 
the specific method for computing the scales, and we also 
provide results obtained by using ground-truth scales, in ad¬ 


dition to the approaches mentioned above. 

EPFL Benchmark. The EPEE Benchmark datasets [27] 
contain from 8 to 30 images, and provide ground-truth ab¬ 
solute motions. 

Results are reported in Tab. 1 and Tab. 2, which show the 
mean errors of motion synchronization before and after ap¬ 
plying a two-step bundle adjustment, as done in [20], where 
in the first step rotations are kept fixed. 

We consider three versions of EIG-SE(3), which dif¬ 
fer for the technique chosen to recover the epipolar scales, 
namely using ground-truth scales (GT), computing scales 
through [3] (MCB), and updating scales iteratively (Iter). 
Our spectral solution is compared with the global SfM 
pipeline described by Moulon et al. [20] and by Ozyesil et 
al. [21]. We also consider the pipeline obtained by combin¬ 
ing the rotation synchronization technique in [4] with the 
translation synchronization method in [7]. As a reference, 
we included in the comparison the sequential SfM pipeline 
Bundler [26]. 

With the exception of Moulon et al. and BUNDLER, for 
which results are taken from [20], all the other methods are 
given the same relative motions as inputs. 

Both EIG-SE(3) and all the analysed techniques achieve 
a high precision, obtaining an average rotation error less 
than 0.1 degrees and an average translation error of the or¬ 
der of millimetres, after the final BA. Our method is able to 
recover camera parameters efficiently, since the motion syn¬ 
chronization step takes about Is for the largest sequences. 

If we concentrate on the EIG-SE(3)-GT columns, we can 
see that it achieves the optimum before BA in most datasets, 
confirming the effectiveness of our method for synchro¬ 
nizing relative motions, when the latter are fully specified. 
Without ground-truth scales, good estimates of motion pa¬ 
rameters are still obtained, and precision increases by using 
MCB rather than the iterative approach. The error after BA 
is always very small and almost equal to the other methods, 
confirming that EIG-SE(3) provides a good starting point 
for bundle adjustment. 

Large-scale Datasets. We test our technique on irregular 
large-scale collections of images taken from [29], for which 
recovering camera orientations/locations is challenging. 

Since our Matlab implementation of Horton’s Algo¬ 
rithm is too slow for large datasets, we do not compute the 
scales through MCB in this experiment. 

We compared EIG-SE(3) with a recent technique - 
called IDSfM [29] - which performs robust translation syn¬ 
chronization. Eollowing the experiments in [29], we used 
the output of Bundler [26] as reference solution, and we 
compute the optimal transformation between this solution 
and our estimate with least median of squares (LMedS), us¬ 
ing correspondences between camera centres. 
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Table 1. Mean angular errors (degrees) on camera rotations for the EPFL benchmark. Moulon et al. is missing in this table because rotation 
errors are not reported in [20]. 


Dataset 

EIG-SE(3)-GT 

EIG-SE(3)-Iter 

EIG-SE(3)-MCB 

Ozyesil et al. 

R-GoDEC-i-Brand et al. 

pre BA 

post BA 

pre BA 

post BA 

pre BA 

post BA 

pre BA 

post BA 

pre BA 

post BA 

HerzJesuPS 

0.04 

0.03 

0.03 

0.03 

0.03 

0.03 

0.06 

0.03 

0.04 

0.03 

HerzJesuP25 

0.06 

0.03 

0.06 

0.04 

0.06 

0.04 

0.14 

0.04 

0.13 

0.04 

FountainPll 

0.03 

0.03 

0.03 

0.04 

0.04 

0.03 

0.03 

0.03 

0.03 

0.03 

EntryPlO 

0.04 

0.02 

0.10 

0.02 

0.11 

0.03 

0.56 

0.04 

0.44 

0.03 

CastleP19 

1.48 

0.06 

1.48 

0.06 

2.46 

0.06 

3.69 

0.05 

1.57 

0.05 

CastlePSO 

0.53 

0.05 

0.47 

0.05 

0.77 

0.05 

1.97 

0.05 

0.78 

0.05 


Table 2. Mean errors (meters) on camera translations for the EPFL benchmark. 


Dataset 

EIG-SE(3)-GT 

EIG-SE(3)-Iter 

EIG-SE(3)-MCB 

Ozyesil et al. 

R-GoDec-i-Brand et al. 

Moulon et al. 

Bundler 

pre BA 

post BA 

pre BA 

post BA 

pre BA 

post BA 

pre BA 

post BA 

pre BA 

post BA 

post BA 

post BA 

HerzJesuPS 

0.004 

0.004 

0.659 

0.004 

0.038 

0.004 

0.007 

0.005 

0.009 

0.004 

0.004 

0.016 

HerzJesuP25 

0.008 

0.008 

1.152 

0.022 

0.357 

0.008 

0.065 

0.009 

0.038 

0.009 

0.005 

0.021 

FountainPll 

0.004 

0.003 

0.236 

0.003 

0.008 

0.003 

0.004 

0.003 

0.006 

0.003 

0.003 

0.007 

EntryPlO 

0.009 

0.008 

0.309 

0.008 

0.349 

0.009 

0.203 

0.010 

0.433 

0.009 

0.006 

0.055 

CastleP19 

0.709 

0.034 

4.986 

0.034 

3.967 

0.035 

1.769 

0.032 

1.493 

0.036 

0.026 

0.344 

CastleP30 

0.212 

0.032 

1.974 

0.035 

3.866 

0.034 

1.393 

0.030 

1.123 

0.030 

0.022 

0.300 


Results are reported in Tab. 3, which shows the median 
errors of motion synchronization before applying bundle 
adjustment. We also report the number of cameras recon¬ 
structed and the percentage of missing pairs, which refer to 
the largest parallel-rigid subgraph, extracted as explained in 
[21]. The results of IDSfM are taken from [29], where rota¬ 
tion errors are not analysed. EIG-SE(3) with iterative scale 
estimate performs equal or better than IDSfM in 7 cases out 
of 11, and it recovers camera rotations accurately. 

Computation times of EIG-SE(3) (Matlab implemen¬ 
tation on a MacBook Air with i5 dual-core @1.3 GHz) re¬ 
ported in Tab. 3 are hardly comparable with those reported 
in [29], as they refer to a compiled code on a much pow¬ 
erful computer. However, if we can assume that the speed 
gain from Matlab to C-H- (for non-trivial algorithms) is at 
least 10 times, as common wisdom suggests, we might then 
conjecture that EIG-SE(3) implemented in C-H- would com¬ 
pare favourably with IDSfM. Moreover, performing paral¬ 
lel computation for updating scales/weights could further 
improve its computational efficiency. 

The rundown of these experiments with real datasets 
shows that, endowed with IRES to withstand outliers and 
combined with a method for estimating the unknown epipo- 
lar scales, EIG-SE(3) can compete with state-of-the-art 
global pipelines. 

4. Conclusion 

We presented a new closed-form method for motion syn¬ 
chronization in SE{3). The method is fast and simple, be¬ 
ing based on a spectral decomposition, and theoretically rel¬ 
evant, for it works in the manifold of rigid motions. 


Table 3. Median errors (rotation in degrees, translation in metres) 
on the datasets from [29] before BA. Boldface denotes the lowest 
translation error. Times are in minutes. 


Dataset 

n 

miss % 

EIG-SE(3)-Iter 

rot. tra. time 

IDSfM 

tra. time 

Roman Forum 

1102 

89 

2.1 

13.5 

16.3 

6.1 

3.5 

Vienna Cathedral 

898 

75 

1.6 

7 

17.6 

6.6 

5 

Alamo 

606 

50 

1.3 

1.5 

14 

1.1 

2.6 

Notre Dame 

553 

32 

0.8 

0.5 

14.7 

10 

2.6 

Tower of London 

489 

81 

2.8 

7.3 

3.0 

11 

1.3 

Montreal N. Dame 

467 

53 

0.6 

O.g 

6.3 

2.5 

1.9 

Yorkminster 

448 

73 

1.9 

7.2 

3.2 

3.4 

2 

Madrid Metropolis 

370 

69 

5 

9.6 

2.5 

9.9 

0.7 

NYC Library 

358 

71 

3.1 

2.5 

2.2 

2.5 

1.3 

Piazza del Popolo 

345 

60 

0.9 

1.6 

2.2 

3.1 

1 

Ellis Island 

240 

33 

0.8 

2.8 

1.8 

3.7 

0.5 


Our experiments showed that our method; i) has the 
same accuracy as its closest competitor [11] but it is much 
faster, and ii) combined with a method for estimating the 
unknown translation norms, it can be profitably used in a 
global structure from motion pipeline with state of the art 
performances. 

The Matlab implementation of E1G-SE(3) will be 
made publicly available. 
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